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We revisit the Hebraud-Lequeux (HL) model for the rheology of jammed materials and argue 
that a possibly important time scale is missing from HL’s initial specification. We show that our 
generalization of the HL model undergoes interesting oscillating instabilities for a wide range of 
parameters, which lead to intermittent, stick-slip flows under constant shear rate. The instability 
we find is akin to the synchronization transition of coupled elements that arises in many different 
contexts (neurons, fireflies, financial bankruptcies, etc.). We hope that our scenario could shed light 
on the commonly observed intermittent, serrated flows of glassy materials under shear. 


I. INTRODUCTION 


Jammed materials are like Tolstoy’s happy families: 
they all behave much in the same way. The (non-linear) 
rheological properties of very different materials such as 
metallic glasses, foams, dense emulsions or granular pack¬ 
ings are remarkably similar. Typically, these materials 
creep very slowly below a certain yield stress, and exhibit 
interesting dynamical features such as “avalanches” [1, 2] 
(perhaps similar, at a different scale, to real earth¬ 
quakes [3]) and intermittent, stick-slip like motion, which 
are still poorly understood - as is the transition from 
“liquid” flows to “jammed” flows [4]. At the microscopic 
scale all of these materials have very different properties 
and can barely be described by one single model. How¬ 
ever on a more coarse grained level, when for example 
considering rheological or plasticity features, these mate¬ 
rials show very similar properties. Such universality has 
motivated researchers to propose simple phenomenologi¬ 
cal models to account for them. It is now well accepted 
that the elementary physical process at play is the yield¬ 
ing of small regions of the material [5-13] (the so-called 
shear transformation zones, STZ), when the local shear 
stress exceeds some (possibly random) threshold. These 
STZ are localized but display long-range elastic interac¬ 
tions [14] that can trigger plastic instabilities. Among 
the most popular approaches to describe the stochastic 
collective evolution of these STZ is the Hebraud-Lequeux 
(HL) model [15], which is a very simple, mean-field de¬ 
scription of the elastic interaction between such STZ, in 
particular how the instability of one of them can trig¬ 
ger the instability of many others (see also Refs. [16-20] 
for alternative descriptions and [21] for a recent review on 
this topic). The HL model, in spite of its simplicity, leads 
to highly non-trivial results, in particular a phase transi¬ 
tion between a jammed, arrested phase and a liquid, flow¬ 
ing phase, with a non-linear Herschel-Bulkley law [22] for 


small shear rate in the jammed phase. This model has 
attracted a renewed interest in the recent years, focusing, 
inter-alia, on a detailed comparison with the competing 
Soft-Glassy-Rheology (SGR) models [23], on avalanche¬ 
like dynamics [24], and on a generalization of HL that 
appropriately accounts for the power-law decay of elastic 
interactions between STZ [25, 26], which deeply modifies 
the mathematical structure of the exact solution of the 
model. 

In this paper, we revisit the original HL model and ar¬ 
gue that some important physical time-scales are missing 
from its original formulation. Although many of the pre¬ 
dictions of HL are unaffected by these modifications, we 
find that an oscillating instability can appear in some re¬ 
gions of parameter space. We find in particular that the 
“liquid” phase at zero shear rate is unstable and becomes, 
close enough the the jammed phase, spontaneously oscil¬ 
lating between a self-sustained liquid and a glass. This in¬ 
stability persists when the shear rate is weak enough, and 
leads, in physical terms, to an intermittent stick-slip flow 
(possibly related to so-called serrated flow). Interest¬ 
ingly, oscillatory instabilities also arise in some variants of 
the SGR models [27], or within simple phenomenological 
constitutive equations for shear-thickening materials [28]. 
However, the underlying physical mechanisms leading to 
instabilities are quite different from the one discussed in 
the present paper. 

The intuition behind our instability in fact comes from 
a very similar model that we recently studied in the 
context of economic crisis “waves” and collective bank 
defaults [29, 30]. In all these cases, the common phe¬ 
nomenon is that the instability/default/bankruptcy of 
a single element can trigger the instabilities of others 
through interactions. The somewhat surprising feature, 
however, in that such a coupling can lead to the syn¬ 
chronization of randomly evolving entities, and genuine 
oscillations - see [30, 31] and references therein. 

The outline of the paper is as follows: we first gener- 
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alize the HL model such as to describe more accurately 
what happens when a single entity fails. We then show 
that the zero shear stationary state computed by HL is 
linearly unstable to some oscillatory mode in a region 
of parameters, in particular close to the jamming transi¬ 
tion. We then extend our calculation to non-zero shear 
rates and find the “phase diagram” of the model, in par¬ 
ticular we highlight the region where stick-slip motion is 
expected. We finally conclude and discuss possible ex¬ 
tensions of our calculations. 


II. THE MODEL 


Let us recall the basic tenets of the HL model [15] and 
introduce the physical items that we believe are missing 
from the original framework. Following HL, the material 
is divided in N regions (or “elements”) labeled by an 
index i = 1, • • • ,N, each characterized by a local stress 
ai- The stress ai of each element performs a random walk 
with: 

• a drift Goj due to the external shear 7 that loads 
elastic energy (Go is the elastic shear modulus), 

• a time-dependent diffusion Dt due to intrinsic noise 
and elastic interactions between elements. 

Whenever \ai\ > the i-th element may become unsta¬ 
ble, and becomes so at a “plastic” rate l/xpi. As soon 
as the element becomes unstable (it is then in a “fluid” 
state), HL postulates that it re-jams (with ct = 0) im¬ 
mediately in a state of zero stress. We rather assume 
that it remains in its fluid state during a typical time 
scale Tfi; in such a state, the element cannot contribute 
to the shear stress (see also Fig. 1 in Ref. [32], where 
Tfl has been called Tei). With rate I/ta, the element 
re-jams with ct = 0, as in HL.^ This is a very impor¬ 
tant ingredient for a correct physical interpretation of 
the model [18, 32]: note that once an element becomes 
unstable, its dynamics cannot be described by the same 
drift-diffusion mechanism. On the contrary, when an ele¬ 
ment becomes unstable, it contributes to the total plastic 
activity at time t, called Ft, and increases the diffusion of 
all other elements identically, in a mean-held way, exactly 
as in the HL model (see [25] for an improved description). 
However, we will assume that the plastic activity at time 
t affects the diffusion coefficient Dt with some delay, in¬ 
stead of instantaneously as HL assumed. This reflects, 
at the mean-held level, the hnite propagation speed of 
information through the sample. 


^ This could actually be generalized to the case where the stress 
drops to a finite fraction of the initial stress (as in [33]), or to a 
non-trivial distribution of initial stresses; the instability reported 
below survives provided the width of that distribution remains 
small. 


A. Mathematical definition of the model 


Our modified model can be described by the same 
Fokker-Planck setting as used by HL [15]. The evolu¬ 
tion of the (unnormalized^) probability P(ct, t) to find an 
element with stress ct at time t is given by: 


P{a,t) = DtP''{a,t) - GojP'{(J,t) 

+ JtS{a) - —P{a,t)H{\a\ - ac) , 

Tpl 

where P[{x) is the Heaviside step function, and with: 
Jt = -{l-<j>t) , 

Tfl 

pOO 


/ oo 

daP{a,t) , 

-00 

Ft = — / da 1 ct 1 P(ct, t) 

J |(7|>cr^ 


( 2 ) 


A = Ant + aw/ ds . 


Here, (pt is the fraction of jammed elements that con¬ 
tribute to the elastic stress and Jt is the flux of elements 
that re-jam between t and t + dt, initially at zero stress. 
One can assume for simplicity^ that the fraction of flu¬ 
idized elements, — contributes a viscous stress pro¬ 
portional to 7 . Hence, the total stress is 

/ OO 

da a P{a, t) + {1 - (j)t)T]j , (3) 

-OO 


where 77 is the microscopic viscosity of the fluid elements. 
Note that satisfies the equation: 

(j)t = —{1 - (j)t) - — [ daP{a,t) . (4) 

^pl J \(7\'> <7 c 

The equation for Dt means that yielding events, when 
they happen, impose an extra random stress onto other 
elements but with some time lag that we model as an 
exponential kernel, with a coupling constant a. is 

an intrinsic noise term (for example temperature) that we 
will choose to be very small in the following. Furthermore 
we choose Gq = 1 without loss of generality^. The con¬ 
trol parameters for this model are fJc, 7 , ta, Tpi, Mnt, 


^ The reason why P((T, t) is not normalised is that it does not 
include the fraction (1 — (pt) of elements that are in the fluid 
state, see Eq. (2). 

^ As it will be clear in the following, the stress St is an observ¬ 
able derived from the model but it does not enter into the main 
equations. Therefore, its precise form does not affect the main 
conclusions of this paper about the phase diagram of the model. 
One could thus choose a different form (e.g. a viscoelastic one) 
for the stress contribution of the fluidized elements. 

Note that Go only appears in Eq. (1) together with 7 . Hence, 
choosing Go = 1 is equivalent to a change of notation, G 07 7 ; 

it is not a choice of units. 
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Because one can fix arbitrarily the units of a and of t, 
there are in fact four plus one independent adimensional 
control parameters: for example 7 := 'yTn/a^ a, Tpi/rfl, 
Q := wTfl, and finally Pe := 'yadDint —> 00 throughout 
this paper. 


Note that in this limit the difference in our definition of P t 
with respect to HL becomes irrelevant, because \a\ = Uc 
is the only contribution to Pj. For simplicity, we also set 
CTc = 1 in the following. 


B. Connection with the Hebraud-Lequeux model 

The standard HL model is recovered if we set I?int = 0, 
Tfl = 0 and w = 00 , i.e: no intrinsic noise, fluid regions 
re-jam instantaneously, and yielding elements affect the 
stress of other regions instantaneously as well. In the 
HL limit rs = 0, we therefore have 0* = 0 and (pt = 1, 
and from Eq. (4) we have Jt = I\a\>a daP{a,t). For 
-Dint = 0 and w —> 00 , we indeed recover HL’s prescrip¬ 
tion Dt = ctPt. Summarizing, we obtain Eq. (1) with 

Jt = — [ daP{a,t) , 

Tpl J\a\>a^ , , 

a f ^ ^ 

Dt = aVt = — / da\a\P{a,t) . 

^pl J\(7\'>(7c 

This coincides with the HL model, with the only differ¬ 
ence that Ft is defined without the \a\. Although this 
does not make any important difference in the proper¬ 
ties of the model, in particular the instability described 
below, we argue that our definition has a more precise 
physical interpretation, since the plastic activity should 
be the average of |cr|. In conclusion, the presence of the 
additional time scales rg and l/w is the main difference 
between our model and HL, and we will see that they 
play a crucial role in the physics of the model. 


C. The limit Tpi 0 

For the ease of analytical calculations, we will set in the 
following Tpi = 0 , which means that an element becomes 
unstable immediately when \a\ > cTc, while keeping a 
finite fluid lifetime Tfl. The HL model corresponds to 
the opposite limit Tpi/rfl —)■ 00 , but it is reasonable to 
expect that in many applications the ratio Tpi/rfl can take 
very different values, from quite small to quite large. In 
this paper we examine analytically the limit Tpi/rfl <C 1 , 
and numerically the regime where Tpi/rfl ^ 1. For small 
Tpi/rfl we find an oscillating instability, which disappears 
when Tpi/rfl exceeds a (parameter dependent) threshold 
value, simply because the synchronisation effect reported 
below cannot set in. So strictly speaking, the instability 
reported here is absent in the original HL setting. 

When Tpi = 0, Eq. (1) is simply complemented by 
an absorbing boundary condition at |(t| = CTc, implying 
P{\(t\ = (Tc) = 0. Furthermore, the plastic activity is 
given by the number of regions that become unstable per 
unit time, hence in Eq. (2) Pj is given by the outgoing 
flux at ±CTc, i.e.: 

Pt = (JcA[D'(-crc,t) - P'dct)] . ( 6 ) 


III. VANISHING SHEAR RATE 

We will consider first the situation with no external 
shear, i.e. 7 = 0 . It turns out that the interesting limit 
corresponds to the case when the intrinsic activity tends 
to zero, i.e. I?int —> 0, as in the HL model. The only 
two adimensional control parameters left are thus a and 

w := WTfl. 

Note that, as we will show in the following (Sect. IV A), 
for Of < 1/2 and = 0 there is a history-dependent 
jammed phase in which the material can self-sustain, in 
the 7 = 0 limit, any finite stress within the interval 
[—(Jc, -t-CTc]. In this case, setting directly 7 = 0 as we do 
below corresponds to a particular choice of history such 
that stationary stress distribution is symmetric and the 
average stress is zero. The main point of this section is 
to illustrate our calculation in the simplest possible case. 


A. The stationary solution: a jamming transition 


We thus look for a stationary solution of Eq. (1) with 
P{a,t) = Poicr), Jt = Jo, etc., which in this case is par¬ 
ticularly simple and reads: 

Po(a > 0) = A(1 - a) , Po{-a) = Pod) , (7) 


with the discontinuity of the slope at the origin related 
to Jo by Eq. (1) as 2DoA = Jq. The fraction of active 
elements <po = fli daPo((j) is simply given by A, so that 
from Eq. (2) 


Jo = -(l-A) ^ A= ^ 

Tfl 2.D0T0 A -1 


( 8 ) 


Eq. ( 6 ) gives Fq = 2ADo and from Eq. (1) we get 
Do = Dint -I- aPo = Pint + 2aADo. We thus obtain a 
self-consistent equation for the stationary value of the 
diffusion coefficient 


2Pq Tfl - 2Po('5a + Ant Tfl) - Pint = 0 , (9) 


where Sa = a — 1/2. This second degree equation has 
two possible solutions that read, in the limit Pint —>■ 0 : 


n- - 
° “ 2<5a 


D+ = 


2ScX CtPintTfl 


2Tfl 


4(5a 


( 10 ) 


where clearly the first solution holds for a < 1/2 and 
the second holds for a > 1/2. When Pint —> 0 there is, 
as first noted by HL, a transition between a fluid phase 
a > 1/2 where activity is self sustained Dq > 0 , and a 
jammed phase a < 1/2 where activity is absent, P/" = 0, 
see Fig. 1 and Ref. [23]. In the jammed phase, one has 
A = 1 and in the liquid phase, A = l/(2a), to leading 
order in Pint- 
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B. Linear stability analysis 

Following [30], we are interested in studying the sta¬ 
bility of the stationary solution with respect to a linear 
perturbation. We thus write P{(J,t) = Po{(t) -\- ePi{(j,t), 
Dt = Dq +eDi{t), <j)t = (j)o +e(f)i{t). Linearizing Eq. (1) 
we get 

Pi{a,t) = DoP{'{a,t)+Di{t)PQ (a)- —(j)i{t)S{a) , (11) 

Tfl 

with the boundary conditions Pi(ct = ±1, t) = 0. There¬ 
fore, 


Eq. (12), taking the derivative with respect to a and 
setting cr = — 1 finally leads to: 



2Do cosh K 


2AV+ —$ 

Tfl 


(18) 


Finally, integrating Eq. (11) over a leads to: 

X<^ = - — ^-2DoIl-2VA. (19) 

Tfl 

In the liquid phase {Dq > 0 and A = l/(2a)) the solv¬ 
ability condition of the system of Eqs. (14), (18), (19) 
relating T>, 11 and $ finally reads, in terms of A = Ata 
and uj = ujTfi' 


pOO pi 

Pi{a,t) = j dr J da G{a,t\d',t — t) 


Di{t - T)Pg{a) - - T)S{a) 

Tfl 


( 12 ) 


where Pl){a) = —2 AS {a), and G is the random walk 
propagator in the interval [— 1 , -fl]. 

Combining Eqs. (2) and ( 6 ), setting Dint = 0 in the 
linear regime we obtain an equation for the diffusion con¬ 
stant Di{t): 


Dtit)=auj [ dsDo[P[{-l,s) - P[{l,s)]e-^‘^*-^^ 

J —OO 

+ 2Aaoj f dsDi{s)e-^d-s) _ 

J —OO 


(13) 


Exponentials of time, with Di{t) = 4’i{t) = ‘he'*'*, 

P{(—l,t) = —P{(l,t) = lie''*, are consistent solutions of 
the linearized equations. We find from Eq. (13): 


V = 


2au}Do ^ 
Cl:( 1 — 2o;^) -f A 


(14) 


In order to compute Pi {a, t) simply, one can notice that 
Ha{a,\) = dTe“^'^G(CT, t|(t, 0) obeys the following 
partial differential equation: 


DoH'd{a, A) - XHg{a, A) = -S{a - a) , (15) 


with boundary condition Hg-{a = ±1,A) = 0. Specializ¬ 
ing to CT = 0 , one finds: 

iL(a>0,A)=l[e-'^"-e-2-e'="] , 

H{-a,X) = H{a,X) 

with = X/Dq. The discontinuity of the first derivative 
at the origin gives: 


2DoA = 


1 

k {\ -I- ’ 


(17) 


so finally H'{—1,X) = —H'{1,X) = l/(2Zlo cosh«:). 
Plugging this in the general expression for Pi{a,t) in 


1 + A = 


w -f A 


w + AcoshC;^ 


or equivalently (for A 7 ^ 0) 


w = 1 — (1 -I- A) cosh \j — . 


( 20 ) 


( 21 ) 


If all the solutions of Eq. (21) have a negative real part, 
the stationary solution is stable, while it is unstable oth¬ 
erwise. This condition defines the transition point add)). 
Note that for a —)■ 00 we obtain from Eq. (21) that 
A = —oj and therefore the stationary state is stable. 


C. Large uj analysis 


Eq. (21) can be analyzed in the limit of large oj (that 
corresponds to the HL case where stress propagation 
is instantaneous throughout the sample), assuming that 
ac{ui) and A are also large (we will see that they are 
proportional to w). We obtain 


w = —A cosh 



( 22 ) 


Let us assume that at the transition there is a single 
unstable mode whose real part crosses zero continuously. 
Thus A = fAi is pure imaginary at Uc- In this case, for 
the r.h.s. of Eq. (22) to be real we need cos i/Xi/Jd<^ = 0 
which implies Ai = a7r^(l -|- 2nd 12 with n G h. Then we 
get 

^ = (-1)” y( 1 + 2n)^ sinh (|(1 + 2n)) . (23) 

Because the system is stable at a —)■ 00 , the stability 
boundary corresponds to the largest value of a for which 
a mode is unstable, and thus to the smallest value of the 
r.h.s. of Eq. (23), which is assumed for n = 0 and gives 


Q:c(^) 

Q 


2 

7r2sinh(f) 


0.088055... , 


(24) 


in excellent agreement with the numerical results of Fig. 1 
when u) is large enough (for example, w = 30 should 
correspond to Oc ~ 2.64). 
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FIG. 1: Top-. Phase diagram of the model in the (0,7) plane 
(recall that 7 = 7rfl) obtained from the solutions of Eq. (44) 
with Di —>■ 0 and different values of 23 = turu. In the liquid 
(LIQ) phase, the system is stationary, and the flow curve is 
basically Newtonian in the whole range of 7. In the jammed 
(JAM) phase, the system has a finite yield stress for 7 —>■ 0. 
Finally, in the oscillating (OSC) phase, the stress is an os¬ 
cillating function of time. Bottom-. Flow curves displaying 
the average stress versus strain for u) = 20, Dint = 0.005, 
and three different values of a corresponding to the jammed, 
oscillating, and liquid phases. Points correspond to numeri¬ 
cal results while lines correspond to analytical results for the 
stationary state. For a = 2.5 the system is liquid and the 
flow curve is Newtonian in all the regimes of 7. For a = 0.25 
the system is jammed for 7 —>■ 0 and the stress tends to a 
finite value (dotted line), but in the simulation due to finite 
Dint we also see a Newtonian regime at small 7 (points and 
dashed line). For a = 1.25 the system is oscillating for small 
7: in this case we report the minimum, average and maximum 
value of the stress (see Fig. 2 for an example). 


D. Summary 

Let us summarize the content of this section. We dis¬ 
cussed the case of vanishing shear rate, 7 = 0 , and intrin¬ 
sic activity, Dint —f 0, assuming Tpi ^ Tfl. The control 
parameters are a and w = wra. We found that: 


• For a > ac(co) « w x 0.088055..., there is a sta¬ 
tionary “liquid” phase with self-sustained activity, 
Dt = Dq > 0 . 


• For 1/2 < a < ac{u}) the liquid phase is in fact 
unstable; in this case, as shown in Fig. 2, the system 
undergoes spontaneous oscillations between periods 
where the number of fluid elements is large and 
periods where it is small. When an external shear is 
applied, this state will correspond to intermittent, 
stick-slip motion. 

Let us emphasize that the liquid phase of the model is 
unphysical in the strict limit Dint = 0 and 7 = 0 , because 
no motion is possible in this case. However, one can think 
that Dint is very small but finite, modeling some external 
noise (e.g. a small temperature), or that 7 —>■ 0+. As we 
discuss next, the more interesting case 7 > 0 is indeed 
perfectly regular in the limit where 7 —>■ 0 . 


IV. NON-ZERO SHEAR RATE 

We now generalize the calculation to a non-zero exter¬ 
nal shear, 7 7 ^ 0, still in the limit Dint 0. In addition 
to a and w = wr^, there is now a third control param¬ 
eter 7 = 7 rfl. We also define ( = -^adDo, which is an 
adimensional “renormalized Peclet number” quantifying 
the importance of drift with respect to the steady-state 
diffusion constant in Eq. (1). Remember however that 
we set CTc = 1 henceforth. 


A. The stationary case 


The stationary solution now reads: 

Po(fT>0)=A(l-e-^(i-'^)) , 

Po{<J < 0) = —Ae“^(l — . 

Integrating Eq. (I) one can relate the discontinuity of 
the slope at the origin to Jo as A = ^ i_|_e-c ■ Moreover, 
the fraction of active elements is (j)o = A (1 — e“^), and 
Jo = (1 — (/o)/'rfl. We thus have three equations for Jo, 
(j)Q and A and the solution is 


7rfi + Jnh(C/2) ’ 

from which we obtain (/q and A thus completing the cal¬ 
culation of Po{cr). 

From Do (o') we obtain the self-consistent equation 
obeyed by Do = 7 /C. Combining Eqs. (2) and (6), we 
have Do = Di„t + oFo = Dint + aDo[D/(-l) - D/(l)]. 
For Dint —>■ 0 we obtain: 

Do 1 a 

7 C 7 + tanh(C/2) 


• For a < 1/2, the system is jammed: no motion 
occurs as Dt = 0 . 


( 27 ) 
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which is more conveniently written as 

7 + tanh(^/ 2 ) = aC, . (28) 


One is now always in a “liquid” phase with a finite dif¬ 
fusion constant Dq > 0 , induced by the shear rate. 

We can also compute the average stress. We obtain 
from Eq. (3) 


_ C-2 + (C + 2)e-^ ^7^ 

+ e~‘^) ^ aC, 


(29) 


where rj = ti/th. Note that for large 7 , Eq. (28) gives 
C ~ 7 /a and thus Eq ~ ^7 is always Newtonian, begin 
dominated by the fluidized elements. A plot of Eq versus 
7 for several values of a is reported in Fig. 1. 

Note that the behavior at small 7 is very different when 
a > 1/2 or when a < 1/2. In the limit 7 —>■ 0: 

• For a > 1/2, it is easy to see graphically that 
Eq. (28) has a unique solution for C for all 7 , and 
that C ~ 7 /( 5 a —)■ 0. Hence, Eq. (25) reduces to 
Eq. (7). Furthermore, 

Do = 7 /(CT-fl) ^ 5a/Tt^ = (30) 

as in Eq. (10), while the stress 

/(si 


is Newtonian, but with a very different prefactor 
from the large 7 regime. 

• For a < 1/2, one finds instead that Eq. (28) has 
multiple solutions when 7 is small enough. Yet at 
large 7 the solution is unique and reducing 7 grad¬ 
ually one sees that V 7 > 0 the largest solution for 
C has to be chosen. This solution for ^ remains fi¬ 
nite for 7 —)• 0. Consequently, Eq. (25) does not 
reduce to Eq. (7), which shows that for a < 1/2 
one cannot set 7 = 0 directly because the limit 
is singular. Setting directly 7 = 0 as in Sec. Ill 
corresponds to choosing a particular solution, cor¬ 
responding to Eq = 0 , which is therefore not the 
physically natural solution, obtained by preparing 
the system under shear and reducing progressively 
7 to zero. In the latter case, Hq oc 7 and Eq re¬ 
mains hnite when 7 —)■ 0 : the system is jammed, 
has no plastic activity and can sustain a finite stress 
even for 7 = 0 . 

Note that C —t 0 for a —)■ 1/2“. Thus Dq = + 0 ( 7 ^) 

for a > 1/2 - see also [23]. 


B. Linear stability analysis 

We linearize once again Eq. (1) with P{(j, t) = Po{cr) + 
ePi{a,t), Dt = Do + eDi{t), (pt = (po + The 


linearized equation reads: 


Pi (a, t) = DoP'^{a, t) - 7 P( (cr, t) 

1 (32) 

+ Di{t)Pi;{a)--P>i{t)5{^) , 

Tfl 

with the boundary conditions Pi{a = ±l,t) = 0. There¬ 
fore 


nOO pi 

Pi{a,t) = j dr j dd — t) 


Di{t - T)Pg{a) - (piit - T)S(a) 

Tfl 


(33) 


where Gj is the biased random walk propagator in the 
interval [—1,-|-1]. The diffusion constant Di{t) satisfies: 

Di{t)=auj [ dsDo[P[i-l,s) - P[il,s)]e-^‘^*-^^ 

J —00 

-I- ^cxuj f (islli(s)e““^*“®^ . 

Do J-00 


(34) 


Assuming exponential (in time) solutions with Di(t) = 
Ve^\ Mt) = T’i'(±l ,0 = Tn±e^‘, we find from 

Eq. (34): 


V = 


aujDo 

w(l — olJo!Do) P a 


(n+ + n_). 


(35) 


Integrating Eq. (32) over a and using Pi{a = ±l,t) = 0 
leads to: 


A$ = -l$-i:)o(n+-fn_)-i?^. (36) 

Tfl Do 

In order to compute Pi (a, t) simply, one can notice that 
Ha{(T,X;j) = dTe~^'^Gj(a,T\d,0) obeys the follow¬ 

ing equation: 

DoH''{a,X;j)-jH'^{a,X;j)-XHg{a,X;'p) = -S(a-d) , 

with boundary condition Hff(i7 = ±1,A;7) = 0. We will 
denote k± the two real roots of the equation: 

DoK^ — — X = 0 , (37) 

and A = K_|_ — AC-. The solution then reads, respectively 
for cr > if and a < d: 

i7s(a,A;7) = A±(a)e"+‘^ + H±(a)e"-" , (38) 

with the following equations coming from boundary con¬ 
ditions: 

(A+ - A_)e'^+'^ -b (5+ - H_)e'"-'^ = 0 , 

k+{A+ - A_)e"+" + Av_(H+ - H_)e"-'^ = . 

Do 


( 39 ) 
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Solving for A± yields: 


A±{a) = 


„=pA-K+cr _ g-K-C 

2Z?o A sinh A 


(40) 


The relevant quantities here are H'^{a = ±1, A; 7 ), which 
are expressed as: 

H'g{a = ±1, A;7) = K+A±{a)e^'^+ + K_B±{a)e^'^- 

= A±((T)Ae^'‘+ ■ 


Using 

DoPoicr) = jPgia) - JoS{a) , 

Pi{a > 0) = P'(a < 0) = , 


in Eq. (33), one finds two contributions. The one coming 
from 6{u) reads: 


= ± 



1 A - e=^'^+ 
— $ - 

Tfl J 2Dq sinh A 


(41) 


while the continuous part gives: 





dcrPo(CT)A±(tT) . 


(42) 


From these we deduce: 


- 




- 





— a - 2.00 


1 1— a = 0.25 



012345678 

t 


FIG. 2: Top: Fraction of jammed elements versns time for a 
population of A = 10 ^ elements undergoing the random walk 
process described by Eq. (1) (see details in the text). Bottom: 
corresponding total stress as defined in Eq. (3). Different 
colors correspond to different values of the coupling constant 
a (black dashed lines correspond to the analytical stationary 
solutions in the stable phase). Other parameters are: 7 = 
0.001, Q = 15, A = 10"® , dt = 10“® and rpi = 0. 


V. NUMERICAL RESULTS 


+ ni 
ui + ui 


JqP , 1 

Do Tfl 

JoVC 


DqK+K- 


) sinh K_ — sinh 

Do sinh A ’ 
(e-«+ - l)(e-«- - 1 ) 
e~^ + 1 


K_l_ sinh K_ — K_ sinh 
sinh A 


(43) 


We now have all the elements to obtain the solvability 
condition of the three equations (35), (36), and (43), re¬ 
lating A, (n_|_ -|-n_) and 4), in terms of A = Ars, w = wta 
and 7 = 7 rfl. We obtain the condition 


5 -C (e-'---l)(e--+-l) ,^ 

wv X + lJ e-'l-bl K+K- 

sinh k_ — k_ sinh 
sinh (k_) — sinh (k+) 

(44) 


The location of the oscillatory instability in the {a, 7 ) 
plane, as predicted by the above equation, forms a kind 
of “bubble” region as shown in Fig. 1. We see that for all 
a S [ 1 / 2 , ac] where the zero-shear state is unstable, there 
exists a critical shear-stress % beyond which the flow is 
stabilized and becomes laminar. Conversely, coming from 
high shear, there is a shear-stress 7 c below which the 
flow becomes intermittent, or stick-slip. This stick-slip 
phenomenon only exists if the elastic coupling between 
elements is neither too large, nor too small. 


The linear stability analysis derived in the previous 
sections is in very good agreement with numerical simu¬ 
lations of the model with a discrete number of elements. 
In order to obtain numerical results we take a set of 
N elements characterized by continuous stress variables 
{a-i(t)}i-i...N and boolean variables {yi{t)}i=i...N (indi¬ 
cating if the element is jammed or is fluidized). At the 
beginning of each time step the ai of jammed elements 
are first independently updated by generating Gaussian 
distributed increments with mean —jdt and round mean 
square ^/2dtDt, where dt is a (sufficiently small) dis¬ 
cretization time step. All elements for which |cri(t)| > Uc 
release their stress with probability dt/rpi, in which case 
they contribute to the Tt+i term in Eq. (2) and become 
fluid. Finally, fluid elements have a probability dt/rfl of 
being re-activated with stress ai = 0. In the limit of large 
N the probability P{a, t) of finding a particle with stress 
between a and a + da evolves according to Eq. (1). 

In Fig. 2 we plot the fraction of jammed elements (j)t 
versus time as well as the total stress Ft for three dif¬ 
ferent values of the coupling constant a with Q = 15, 
Tpi = 0, 7 = 0.001. One clearly observes that the liquid 
and jammed phases are separated by a stick-slip regime 
for intermediate values of the coupling constant a, as 
predicted analytically. A more extensive numerical ex¬ 
ploration of the model (not shown here) shows a very 
good agreement with analytical results for the location 
of phase boundaries depicted in Fig. 1. 

In order to have a better understanding of the dynam- 





























FIG. 3: Empirical distribution of stress variables ct; for a pop¬ 
ulation of = 10® elements with zero shear rate. Different 
lines are relative to different snapshots taken at equally spaced 
time intervals of duration 0.01 (the corresponding value for 
the fraction of jammed elements as a function of time is plot¬ 
ted in the inset). Red lines correspond to the build up of the 
instability (the lower line is at t = 0.03 for example, approx¬ 
imately corresponding to the minimum of 0(t) in the inset) 
while black lines correspond to elements progressively relax¬ 
ing at (T = 0. For illustrative purposes we set here 7 = 0 
while other parameters are: a = 1 , th = 0 .01, Di = 0.01, 
dt = 10“®, tD = 30 and Tpi = 0. 

ics of the model in the unstable regime it is instructive 
to look at the empirical probability distribution P((t, t) 
at different times corresponding to the accelerated build¬ 
up of stress and its subsequent relaxation. This is done 
in Fig. 3 where we plot the empirical P(a,t) at snap¬ 
shots equally spaced in time (see the caption for details). 
We observe that initially, P{a,t) is peaked around zero. 
Then, as some elements become unstable, the stress dif¬ 
fusion constant increases, bringing more elements close 
to (Tc- This feedback can become explosive and lead to a 
“spike” in the number of elements that become unstable 
simultaneously, see Fig. 3, inset. 

Finally, since all the calculations above where made 
assuming Tpi = 0, we have explored numerically how the 
instability behaves when Tpi > 0. In this case, we find 
that the qualitative features of Fig. 1 are left unchanged 
as long as Tpi < where the value of depends on 
other parameters, whereas the instability disappears al¬ 
together when Tpi > Tpj. For example, in the zero shear 
rate case, we End for w = 30 and a = 2 that the insta¬ 
bility disappears at r*j « 0.857 rg. The reason is quite 
simple: since jammed elements become unstable at ran¬ 
dom times with a dispersion ^ Tpi, a large Tpi maims the 
synchronisation mechanism that leads to the oscillating 


instability discovered here. Since the HL model corre¬ 
sponds to Tpi/rfl —>■ oo, this oscillating instability would 
not have shown up in the original HL setting. Note that 
similarly, a small w (corresponding to a large dispersion 
in the arrival time of the stress signal) is preventing the 
instability to occur. 

VI. CONCLUSION 

In this paper, we have revisited the now classic 
Hebraud-Lequeux model for the rheology of jammed ma¬ 
terials. We have argued that a possibly important time 
scale is missing from HL’s initial specification, namely the 
characteristic time for a fluidized element to re-jam (cho¬ 
sen to be zero in HL), whose importance has also been 
stressed in Refs. [18, 20, 32]. A linear stability analysis of 
the generalized HL model shows that the steady-state so¬ 
lution is in fact unstable for a wide range of parameters, 
and leads to intermittent, stick-slip flows under a con¬ 
stant shear rate. The stick-slip motion disappears when 
the shear-rate is large (as expected), or when the time 
for a jammed element to become unstable is large com¬ 
pared to the re-jamming time, or else when the elastic 
coupling between jammed elements is either too weak or 
too strong (see Fig. 1 above). 

The instability we End is akin to the synchronization 
transition of coupled elements that arises in many dif¬ 
ferent contexts (neurons, EreEies, Enancial bankruptcies, 
etc.), see [30, 31] and references therein. Similar instabil¬ 
ities are found in some variants of the SGR models [27] 
and within simple phenomenological constitutive equa¬ 
tions for shear-thickening materials [28], although the 
underlying physical mechanism does not seem to be re¬ 
lated to the synchronisation eEect discussed here. We 
hope that our scenario could shed light on the commonly 
observed intermittent, serrated flows of glassy materials 
under shear. It would also be quite interesting to study 
our instability in finite dimensions (rather than in the HL 
mean-held limit), and along the lines recently suggested 
by Lin and Wyart [25], where the HL diffusion term is 
replaced by a long-range. Levy flight diflusion. It is not 
immediately clear whether the stick-slip instability sur¬ 
vives such a deep modiflcation of the mathematical and 
physical nature of the model. 
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